Spatial Ecology and Macroecology

Practical - Week 3

Florencia Grattarola
grattarola@fzp.czu.cz

MOBI Lab

2024-10-14

What are we going to see today?

We will explore the potential drivers of species distribution.

  1. Data sources (climate, habitat types, land cover, socioeconomic)
  2. Data download and processing
  3. Visualisation/mapping
  4. Correlations between predictors
  5. Description of patterns (vs. biodiversity)

Exercises:

  1. Peru: drivers of amphibians richness patterns 🐸
  2. France: drivers of trees richness patterns 🌳

Data sources

worldclim.org

WorldClim Is a database of high spatial resolution global weather and climate data. You can download gridded weather and climate data for historical (near current) and future conditions.

chelsa-climate.org

CHELSA (Climatologies at high resolution for the earth’s land surface areas) is a very high resolution (30 arc sec, ~1km) global downscaled climate data set. It includes climate layers for various time periods and variables, ranging from the Last Glacial Maximum, to the present, to several future scenarios.

nature.com/articles/s41597-020-00599-8

A global, spatially explicit characterization of 47 terrestrial habitat types, as defined in the International Union for Conservation of Nature (IUCN) habitat classification scheme.

Habitatmapping: https://github.com/Martin-Jung/Habitatmapping

ecoregions.appspot.com

The RESOLVE Ecoregions dataset, updated in 2017, offers a depiction of the 846 terrestrial ecoregions that represent our living planet.

land.copernicus.eu/pan-european/corine-land-cover

The CORINE Land Cover (CLC) consists of an inventory of land cover in 44 classes (from Europe). The inventory was initiated in 1985. Updates have been produced in 2000, 2006, 2012, and 2018.

modis.gsfc.nasa.gov

Terra MODIS and Aqua MODIS satellites are viewing the entire Earth’s surface every 1 to 2 days, acquiring data in 36 spectral bands or groups of wavelengths.

MODIStsp: https://github.com/ropensci/MODIStsp/

sedac.ciesin.columbia.edu

SEDAC, the Socioeconomic Data and Applications Center, focuses on human interactions in the environment and seeks to develop and operate applications that support the integration of socioeconomic and earth science data.

EXCERCISE Peru: drivers of amphibians richness

Peru: drivers of amphibians richness 🐸

Steps:

  1. Get a map of Peru
  2. Get 5 different drivers (raster)
  3. Process all the layers to the same resolution and the extent of Peru
  4. Visualise predictors
  5. Calculate the correlation between the different predictors (and compare with biodiversity)

Some preparation before starting to code

Biodiversiy data download

Please download the following file and store it on you data folder:

  • amphibians_grid.rds

To create this file, I followed steps from the previous practicals. First I downloaded and cleaned data from GBIF, and then I created a grid cell for Peru of 50 x 50 km.

Check download_amphibiansPE_data_from_GBIF.R to see how the data were downloaded, and preparing_amphibiansPE_grid.R to see how I created the grid

Biodiversiy data download

This is how the data look like

Packages

We will use sf, rnaturalearth, and tidyverse.

library(pacman)
pacman::p_load(sf,  # Manage spatial data as simple features dataframes
               rnaturalearth, # useful base maps
               tidyverse) # install the package if you haven't already

sf_use_s2(FALSE) # switch spherical geometry off

Packages

We will also use terra to manage and analyse rasters and vectors, tmap for plotting maps, geodata to get multiple sources of spatial data, and MODIStsp for processing time series from MODIS data.

pacman::p_load(terra, # Work with rasters and vectorial spatial data
               tmap, # Cool maps
               GGally, # Cool correlation maps
               geodata, # Get data from multiple sources
               MODIStsp) # Get data from MODIS

Peru: drivers of amphibians richness 🐸

  1. Get a map of Peru

We will need,

  • The polygon and spatial extent of Peru
  • A planar projection of the area.

Peru: drivers of amphibians richness 🐸

  1. Get a map of Peru
peru <- ne_countries(
  country = "Peru",
  returnclass = "sf" # We want an object of class sf
)

Let’s see how it looks.

Peru: drivers of amphibians richness 🐸

  1. Get a map of Peru
ggplot() + 
  geom_sf(data=peru, fill='white') +
  coord_sf()

Peru: drivers of amphibians richness 🐸

  1. Get a map of Peru

To get a planar projection of the area (visit https://projectionwizard.org/).

proj_peru <- 'PROJCS["ProjWiz_Custom_Transverse_Cylindrical_Equal_Area",
 GEOGCS["GCS_WGS_1984",
  DATUM["D_WGS_1984",
   SPHEROID["WGS_1984",6378137.0,298.257223563]],
  PRIMEM["Greenwich",0.0],
  UNIT["Degree",0.0174532925199433]],
 PROJECTION["Transverse_Cylindrical_Equal_Area"],
 PARAMETER["False_Easting",0.0],
 PARAMETER["False_Northing",0.0],
 PARAMETER["Central_Meridian",-74.7949219],
 PARAMETER["Scale_Factor",1.0],
 PARAMETER["Latitude_Of_Origin",0.0],
 UNIT["Meter",1.0]]'

We will use it later.

Peru: drivers of amphibians richness 🐸

  1. Get 5 different drivers (raster)
  • Annual mean temperature (WorldClim BIO1)
  • Precipitation seasonality (WorldClim BIO15)
  • Elevation (SRTM)
  • Land-cover (MODIS Terra MCD12Q1)
  • Human footprint index (SEDAC)

Peru: drivers of amphibians richness 🐸

  1. Get 5 different drivers (raster)

When working with rasters, there are three basic procedures we can recommend you:

  1. Crop: Cut the rasters to the geographic extent you need. This is a square cut.
  2. Reproject: Make sure that all your rasters are on the same projection. But be aware that by reprojecting, you are introducing an additional source of error.
  3. Mask: Set to NA all the cells outside your polygon(s) of interest.

Peru: drivers of amphibians richness 🐸

  1. Get 5 different drivers (raster)

We will create a variable with the spatial extent of Peru to use it to crop the raster data.

round(st_bbox(peru),0) # Check the extent of the country (round values)
xmin ymin xmax ymax 
 -81  -18  -69    0 
round(st_bbox(peru),0) + c(-1,-1,1,1) #This way we add space around the country
xmin ymin xmax ymax 
 -82  -19  -68    1


Now let’s create an object with the extent.

extent <- round(st_bbox(peru),0) + c(-1,-1,1,1) # Bounding box plus 1 degree

Peru: drivers of amphibians richness 🐸

  1. Get 5 different drivers (raster)

We will also set the path were we are going to save data.

geodata_path("data/")


This is required to use the package geodata. From now on, every time you try to get a dataset with geodata functions, the package will check if the data already exists in your data folder and load it from there.

Peru: drivers of amphibians richness 🐸

  1. Get 5 different drivers (raster): climate

Let’s first download the WorldClim data (version 2.1 climate data for 1970-2000):

  • Annual mean temperature (bio1)
  • Annual precipitation (bio15)


You can download these data from https://www.worldclim.org/data/worldclim21.html
or use the package geodata.

Peru: drivers of amphibians richness 🐸

  1. Get 5 different drivers (raster): climate

We will use the function worldclim_country() from geodata.

tavg <- worldclim_country(var = "tavg",country = "PER") # PER is the three letters ISO code for PERU
tavg <- crop(tavg, extent)
tavg
class       : SpatRaster 
dimensions  : 2220, 1560, 12  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : -81.5, -68.5, -18.5, 0  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : PER_wc2.1_30s_tavg.tif 
names       : PER_w~avg_1, PER_w~avg_2, PER_w~avg_3, PER_w~avg_4, PER_w~avg_5, PER_w~avg_6, ... 
min values  :        -5.4,        -5.4,        -5.2,        -5.8,        -6.9,        -8.6, ... 
max values  :        27.4,        27.3,        27.3,        27.3,        27.0,        26.8, ...

Peru: drivers of amphibians richness 🐸

  1. Get 5 different drivers (raster): climate

Now get the precipitation data from Worldclim.

prec <- worldclim_country(
  var = "prec",
  country = "PER"
)

prec <- crop(prec, extent)

prec
class       : SpatRaster 
dimensions  : 2220, 1560, 12  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : -81.5, -68.5, -18.5, 0  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : PER_wc2.1_30s_prec.tif 
names       : PER_w~rec_1, PER_w~rec_2, PER_w~rec_3, PER_w~rec_4, PER_w~rec_5, PER_w~rec_6, ... 
min values  :           0,           0,           0,           0,           0,           0, ... 
max values  :         973,         832,         617,         499,         477,         516, ...

Peru: drivers of amphibians richness 🐸

  1. Get 5 different drivers (raster): elevation

We will use the function elevation_30s() from geodata. These data come from the Shuttle Radar Topographic Mission and are aggregated to 1km resolution.

elev <- elevation_30s(country = "PER")

elev <- crop(elev, extent)

elev
class       : SpatRaster 
dimensions  : 2232, 1560, 1  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : -81.5, -68.5, -18.5, 0.1  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : PER_elv_msk.tif 
name        : PER_elv_msk 
min value   :         -25 
max value   :        6547

Peru: drivers of amphibians richness 🐸

  1. Get 5 different drivers (raster): land-cover

Now, let’s get the landcover data. For this we will use the package MODIStsp.

MODIStsp has not been maintained

The package may not be working properly.

MODIStsp_get_prodnames()
  [1] "Surf_Ref_8Days_500m (M*D09A1)"                   
  [2] "Surf_Ref_Daily_005dg (M*D09CMG)"                 
  [3] "Surf_Ref_Daily_500m (M*D09GA)"                   
  [4] "Surf_Ref_Daily_250m (M*D09GQ)"                   
  [5] "Surf_Ref_8Days_250m (M*D09Q1)"                   
  [6] "Snow_Cov_Daily_500m (M*D10A1)"                   
  [7] "Snow_Cov_8-Day_500m (M*D10_A2)"                  
  [8] "Snow_Cov_Day_0.05Deg (M*D10C1)"                  
  [9] "Snow_Cov_8-Day0.05Deg CMG (M*D10C2)"             
 [10] "Snow_Cov_Month_0.05Deg CMG (M*D10CM)"            
 [11] "Surf_Temp_Daily_005dg (M*D11C1)"                 
 [12] "Surf_Temp_Daily_1Km (M*D11A1)"                   
 [13] "Surf_Temp_8Days_1Km (M*D11A2)"                   
 [14] "Surf_Temp_Daily_GridSin (M*D11B1)"               
 [15] "Surf_Temp_8Days_GridSin (M*D11B2)"               
 [16] "Surf_Temp_Monthly_GridSin (M*D11B3)"             
 [17] "Surf_Temp_8Days_005dg (M*D11C2)"                 
 [18] "Surf_Temp_Monthly_005dg (M*D11C3)"               
 [19] "LST_3band_emissivity_Daily_1km (M*D21A1D)"       
 [20] "LST_3band_emissivity_Daily_1km_night (M*D21A1N)" 
 [21] "LST_3band_emissivity_8day_1km (M*D21A2)"         
 [22] "BRDF_Albedo_ModelPar_Daily_500m (MCD43A1)"       
 [23] "BRDF_Albedo_Quality_Daily_500m (MCD43A2)"        
 [24] "Albedo_Daily_500m (MCD43A3)"                     
 [25] "BRDF_Adj_Refl_Daily_500m (MCD43A4)"              
 [26] "BRDF_Albedo_ModelPar_Daily_005dg (MCD43C1)"      
 [27] "BRDF_Albedo_Quality_Daily_005dg (MCD43C2)"       
 [28] "Albedo_Daily_005dg (MCD43C3)"                    
 [29] "BRDF_Adj_Refl_16Day_005dg (MCD43C4)"             
 [30] "AlbPar_1_B1_Daily_30ArcSec (MCD43D01)"           
 [31] "AlbPar_2_B1_Daily_30ArcSec (MCD43D02)"           
 [32] "AlbPar_3_B1_Daily_30ArcSec (MCD43D03)"           
 [33] "AlbPar_1_B2_Daily_30ArcSec (MCD43D04)"           
 [34] "AlbPar_2_B2_Daily_30ArcSec (MCD43D05)"           
 [35] "AlbPar_3_B2_Daily_30ArcSec (MCD43D06)"           
 [36] "AlbPar_1_B3_Daily_30ArcSec (MCD43D07)"           
 [37] "AlbPar_2_B3_Daily_30ArcSec (MCD43D08)"           
 [38] "AlbPar_3_B3_Daily_30ArcSec (MCD43D09)"           
 [39] "AlbPar_1_B4_Daily_30ArcSec (MCD43D10)"           
 [40] "AlbPar_2_B4_Daily_30ArcSec (MCD43D11)"           
 [41] "AlbPar_3_B4_Daily_30ArcSec (MCD43D12)"           
 [42] "AlbPar_1_B4_Daily_30ArcSec (MCD43D13)"           
 [43] "AlbPar_2_B4_Daily_30ArcSec (MCD43D14)"           
 [44] "AlbPar_3_B4_Daily_30ArcSec (MCD43D15)"           
 [45] "AlbPar_1_B5_Daily_30ArcSec (MCD43D16)"           
 [46] "AlbPar_2_B5_Daily_30ArcSec (MCD43D17)"           
 [47] "AlbPar_3_B5_Daily_30ArcSec (MCD43D18)"           
 [48] "AlbPar_1_B6_Daily_30ArcSec (MCD43D19)"           
 [49] "AlbPar_2_B6_Daily_30ArcSec (MCD43D20)"           
 [50] "AlbPar_3_B6_Daily_30ArcSec (MCD43D21)"           
 [51] "AlbPar_1_Vis_Daily_30ArcSec (MCD43D22)"          
 [52] "AlbPar_2_Vis_Daily_30ArcSec (MCD43D23)"          
 [53] "AlbPar_3_Vis_Daily_30ArcSec (MCD43D24)"          
 [54] "AlbPar_1_NIR_Daily_30ArcSec (MCD43D25)"          
 [55] "AlbPar_2_NIR_Daily_30ArcSec (MCD43D26)"          
 [56] "AlbPar_3_NIR_Daily_30ArcSec (MCD43D27)"          
 [57] "AlbPar_1_SWIR_Daily_30ArcSec (MCD43D28)"         
 [58] "AlbPar_2_SWIR_Daily_30ArcSec (MCD43D29)"         
 [59] "AlbPar_3_SWIR_Daily_30ArcSec (MCD43D30)"         
 [60] "BRDF_Albedo_Quality_Daily_30ArcSec (MCD43D31)"   
 [61] "BRDF_Albedo_SolNoon_Daily_30ArcSec (MCD43D32)"   
 [62] "Alb_ValObs_B1_Daily_30ArcSec (MCD43D33)"         
 [63] "Alb_ValObs_B2_Daily_30ArcSec (MCD43D34)"         
 [64] "Alb_ValObs_B3_Daily_30ArcSec (MCD43D35)"         
 [65] "Alb_ValObs_B4_Daily_30ArcSec (MCD43D36)"         
 [66] "Alb_ValObs_B5_Daily_30ArcSec (MCD43D37)"         
 [67] "Alb_ValObs_B6_Daily_30ArcSec (MCD43D38)"         
 [68] "Alb_ValObs_B7_Daily_30ArcSec (MCD43D39)"         
 [69] "BRDF_Albedo_Snow_Daily_30ArcSec (MCD43D40)"      
 [70] "BRDF_Alb_Unc_Daily_30ArcSec (MCD43D41)"          
 [71] "BRDF_Alb_BSA_B1_Daily_30ArcSec (MCD43D42)"       
 [72] "BRDF_Alb_BSA_B2_Daily_30ArcSec (MCD43D43)"       
 [73] "BRDF_Alb_BSA_B3_Daily_30ArcSec (MCD43D44)"       
 [74] "BRDF_Alb_BSA_B4_Daily_30ArcSec (MCD43D45)"       
 [75] "BRDF_Alb_BSA_B5_Daily_30ArcSec (MCD43D46)"       
 [76] "BRDF_Alb_BSA_B6_Daily_30ArcSec (MCD43D47)"       
 [77] "BRDF_Alb_BSA_B7_Daily_30ArcSec (MCD43D48)"       
 [78] "BRDF_Alb_BSA_Vis_Daily_30ArcSec (MCD43D49)"      
 [79] "BRDF_Alb_BSA_NIR_Daily_30ArcSec (MCD43D50)"      
 [80] "BRDF_Alb_BSA_SWIR_Daily_30ArcSec (MCD43D51)"     
 [81] "BRDF_Alb_WSA_B1_Daily_30ArcSec (MCD43D52)"       
 [82] "BRDF_Alb_WSA_B2_Daily_30ArcSec (MCD43D53)"       
 [83] "BRDF_Alb_WSA_B3_Daily_30ArcSec (MCD43D54)"       
 [84] "BRDF_Alb_WSA_B4_Daily_30ArcSec (MCD43D55)"       
 [85] "BRDF_Alb_WSA_B5_Daily_30ArcSec (MCD43D56)"       
 [86] "BRDF_Alb_WSA_B6_Daily_30ArcSec (MCD43D57)"       
 [87] "BRDF_Alb_WSA_B7_Daily_30ArcSec (MCD43D58)"       
 [88] "BRDF_Alb_WSA_Vis_Daily_30ArcSec (MCD43D59)"      
 [89] "BRDF_Alb_WSA_NIR_Daily_30ArcSec (MCD43D60)"      
 [90] "BRDF_Alb_WSA_SWIR_Daily_30ArcSec (MCD43D61)"     
 [91] "BRDF_Albedo_NBAR_Band1_Daily_30ArcSec (MCD43D62)"
 [92] "BRDF_Albedo_NBAR_Band2_Daily_30ArcSec (MCD43D63)"
 [93] "BRDF_Albedo_NBAR_Band3_Daily_30ArcSec (MCD43D64)"
 [94] "BRDF_Albedo_NBAR_Band4_Daily_30ArcSec (MCD43D65)"
 [95] "BRDF_Albedo_NBAR_Band5_Daily_30ArcSec (MCD43D66)"
 [96] "BRDF_Albedo_NBAR_Band6_Daily_30ArcSec (MCD43D67)"
 [97] "BRDF_Albedo_NBAR_Band7_Daily_30ArcSec (MCD43D68)"
 [98] "Vegetation_Indexes_16Days_500m (M*D13A1)"        
 [99] "Vegetation_Indexes_16Days_1Km (M*D13A2)"         
[100] "Vegetation_Indexes_Monthly_1Km (M*D13A3)"        
[101] "Vegetation_Indexes_16Days_005dg (M*D13C1)"       
[102] "Vegetation_Indexes_Monthly_005dg (M*D13C2)"      
[103] "Vegetation Indexes_16Days_250m (M*D13Q1)"        
[104] "LAI_8Days_500m (MCD15A2H)"                       
[105] "LAI_4Days_500m (MCD15A3H)"                       
[106] "LAI_8Days_500m (M*D15A2H)"                       
[107] "Net_ET_8Day_500m (M*D16A2)"                      
[108] "Net_ETgf_8Day_500m (M*D16A2GF)"                  
[109] "Net_ETgf_Yearly_500m (M*D16A3GF)"                
[110] "Gross_PP_8Days_500m (M*D17A2H)"                  
[111] "Gross_PP_GapFil_8Days_500m (M*D17A2HGF)"         
[112] "Net_PP_GapFil_Yearly_500m (M*D17A3HGF)"          
[113] "Veg_Cont_Fields_Yearly_250m (MOD44B)"            
[114] "Land_Wat_Mask_Yearly_250m (MOD44W)"              
[115] "Burned_Monthly_500m (MCD64A1)"                   
[116] "ThermalAn_Fire_Daily_1Km (M*D14A1)"              
[117] "ThermalAn_Fire_8Days_1Km (M*D14A2)"              
[118] "LandCover_Type_Yearly_005dg (MCD12C1)"           
[119] "LandCover_Type_Yearly_500m (MCD12Q1)"            
[120] "LandCover_Dynamics_Yearly_500m (MCD12Q2)"        
[121] "Dwnwrd_Srw_Rad_3h_005dg (MCD18A1)"               
[122] "Dwnwrd_PAR_3h_005dg (MCD18A2)"                   
[123] "MAIA_Land_Surf_BRF (MCD19A1)"                    
[124] "MAIA_Land_AOT_daily (MCD19A2)"

Peru: drivers of amphibians richness 🐸

  1. Get 5 different drivers (raster): land-cover

To download MODIS data through the NASA http server, we need to create a profile at https://urs.earthdata.nasa.gov/home to get a user and password.


Today, I’ll present you an example and provide you with the data already processed, but you can do this at home.


We will use the Land Cover Type Yearly L3 Global 500m, and download data for Peru from the year 2020.

Peru: drivers of amphibians richness 🐸

  1. Get 5 different drivers (raster): land-cover

With MODIStsp_get_prodlayers() you can see all the layers of the product:

MODIStsp_get_prodlayers("MCD12Q1", prodver = "061")$bandnames
 [1] "LC1"                 "LC2"                 "LC3"                
 [4] "LC4"                 "LC5"                 "LC_Prop1"           
 [7] "LC_Prop2"            "LC_Prop3"            "LC_Prop1_Assessment"
[10] "LC_Prop2_Assessment" "LC_Prop3_Assessment" "LC_QC"              
[13] "LC_LW"


There are five different land cover classification schemes, we will be using the primary land cover scheme (LC1) which identifies 17 classes defined by the IGBP (International Geosphere-Biosphere Programme), including 11 natural vegetation classes, 3 human-altered classes, and 3 non-vegetated classes.

Peru: drivers of amphibians richness 🐸

  1. Get 5 different drivers (raster): land-cover

We will use the function MODIStsp() to download data:

MODIStsp(
  ...,
  gui = TRUE,
  out_folder = NULL,
  out_folder_mod = NULL,
  opts_file = NULL,
  selprod = NULL,
  prod_version = NULL,
  bandsel = NULL,
  quality_bandsel = NULL,
  indexes_bandsel = NULL,
  sensor = NULL,
  download_server = NULL,
  downloader = NULL,
  user = NULL,
  password = NULL,
  download_range = NULL,
  start_date = NULL,
  end_date = NULL,
  spatmeth = NULL,
  start_x = NULL,
  end_x = NULL,
  start_y = NULL,
  end_y = NULL,
  bbox = NULL,
  spafile = NULL,
  out_projsel = NULL,
  output_proj = NULL,
  out_res_sel = NULL,
  out_res = NULL,
  resampling = NULL,
  reprocess = NULL,
  delete_hdf = NULL,
  nodata_change = NULL,
  scale_val = NULL,
  ts_format = NULL,
  out_format = NULL,
  compress = NULL,
  test = NULL,
  n_retries = 5,
  verbose = TRUE,
  parallel = TRUE
)

Peru: drivers of amphibians richness 🐸

  1. Get 5 different drivers (raster): land-cover
MODIStsp(
  gui = FALSE,
  out_folder = "data/",
  out_folder_mod = "data/",
  selprod = "LandCover_Type_Yearly_500m (MCD12Q1)",
  prod_version = "061",
  bandsel = "LC1",
  sensor = "Terra",
  user = your_user, # your username for NASA http server
  password = your_pass, # your password for NASA http server
  start_date = "2020.01.01",
  end_date = "2020.12.31",
  verbose = TRUE,
  bbox = extent, # bbox covering Peru
  spatmeth = "bbox",
  out_format = "GTiff",
  compress = "LZW",
  out_projsel = "User Defined",
  output_proj = 4326,
  delete_hdf = TRUE,
  parallel = TRUE
)

Peru: drivers of amphibians richness 🐸

  1. Get 5 different drivers (raster): land-cover

Let’s read the layer and process it.

land <- rast("data/MCD12Q1.061_LC_Type1_doy2020001_aid0001.tif")

land <- crop(land , extent)

land
class       : SpatRaster 
dimensions  : 4800, 3360, 1  (nrow, ncol, nlyr)
resolution  : 0.004166667, 0.004166667  (x, y)
extent      : -82, -68, -19, 1  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source(s)   : memory
varname     : MCD12Q1.061_LC_Type1_doy2020001_aid0001 
name        : MCD12Q1.061_LC_Type1_doy2020001_aid0001 
min value   :                                       1 
max value   :                                      17

Peru: drivers of amphibians richness 🐸

  1. Get 5 different drivers (raster): anthropogenic

Again using the package geodata.

footp <- footprint(year = 2009)

footp <- crop(footp, extent)

footp
class       : SpatRaster 
dimensions  : 2400, 1680, 1  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : -82, -68, -19, 1  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source(s)   : memory
varname     : wildareas-v3-2009-human-footprint_geo 
name        : wildareas-v3-2009-human-footprint 
min value   :                                 0 
max value   :                                46

Peru: drivers of amphibians richness 🐸

  1. Process all the layers to the same resolution and the extent of Peru

Rasters should have the same spatial extent and resolution… If they don’t match… what should we do?

compareGeom(tavg, prec, elev, footp, land, stopOnError = F, messages = T, res = T)
[1] FALSE

Peru: drivers of amphibians richness 🐸

  1. Process all the layers to the same resolution and the extent of Peru
st_bbox(tavg)
 xmin  ymin  xmax  ymax 
-81.5 -18.5 -68.5   0.0 
st_bbox(prec) 
 xmin  ymin  xmax  ymax 
-81.5 -18.5 -68.5   0.0 
st_bbox(elev) 
 xmin  ymin  xmax  ymax 
-81.5 -18.5 -68.5   0.1 
st_bbox(footp)
xmin ymin xmax ymax 
 -82  -19  -68    1 
st_bbox(land)
xmin ymin xmax ymax 
 -82  -19  -68    1 
#Compare with the country extent
st_bbox(peru)
       xmin        ymin        xmax        ymax 
-81.4109426 -18.3479754 -68.6650797  -0.0572055

Peru: drivers of amphibians richness 🐸

  1. Process all the layers to the same resolution and the extent of Peru

We will use the package terra to process the rasters.
First we need to create a new raster to use as a template.

Bug in terra

We could do:

r <- rast(res = 10000,
          ext = extention, # or even `ext(vect(peru))`
          crs = "epsg:4326")

But there’s a bug :O
I reported it https://github.com/rspatial/terra/issues/1618

Peru: drivers of amphibians richness 🐸

  1. Process all the layers to the same resolution and the extent of Peru

We will use the package terra to process the rasters.
First we need to create a new raster to use as a template.

r <- rast(res = 10000, # Target cell size
          ext(vect(st_transform(peru, proj_peru))), # Extent of the raster
          crs = proj_peru) %>% # Useful planar projection
  project(., "epsg:4326") # Projection of our rasters

r
class       : SpatRaster 
dimensions  : 203, 147, 1  (nrow, ncol, nlyr)
resolution  : 0.09028797, 0.09028797  (x, y)
extent      : -81.74007, -68.46773, -18.40142, -0.07296261  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)

Peru: drivers of amphibians richness 🐸

  1. Process all the layers to the same resolution and the extent of Peru

We will resample() all of the rasters to our target resolution.

tavg <- resample(tavg, r, method = "bilinear") # for continuous data!
prec <- resample(prec, r, method = "bilinear")
elev <- resample(elev, r, method = "bilinear")
footp <- resample(footp, r, method = "near")
land <- resample(land, r, method = "near")  # for categorical data!

And we can compare the rasters after resampling.

compareGeom(tavg, prec, elev, footp, land, stopOnError = F, messages = T, res = T)
[1] TRUE

Peru: drivers of amphibians richness 🐸

  1. Process all the layers to the same resolution and the extent of Peru

And finally we will use mask() to intersect the layer with Peru.

tavg <- mask(tavg, vect(peru))
prec <- mask(prec, vect(peru))
elev <- mask(elev, vect(peru))
footp <- mask(footp, vect(peru))
land <- mask(land, vect(peru))

Peru: drivers of amphibians richness 🐸

  1. Visualise predictors

We will use the package tmap.

tmap_mode(mode = "view")

Peru: drivers of amphibians richness 🐸

  1. Visualise predictors

Annual mean temperature (WorldClim BIO1)

Peru: drivers of amphibians richness 🐸

  1. Visualise predictors

Precipitation seasonality (WorldClim BIO15)

Peru: drivers of amphibians richness 🐸

  1. Visualise predictors

Elevation (WorldClim)

Peru: drivers of amphibians richness 🐸

  1. Visualise predictors

Land-cover (MODIS Terra MCD12Q1). We need to make sure values are factors.

names(land)
[1] "MCD12Q1.061_LC_Type1_doy2020001_aid0001"
land <- as.factor(land$MCD12Q1.061_LC_Type1_doy2020001_aid0001)

Peru: drivers of amphibians richness 🐸

  1. Visualise predictors

Then we can add values to the levels.

levels(land[[1]]) <- data.frame(
  ID = c(1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 17),
  label = c(
    "Evergreen Needleleaf Forest",
    "Evergreen Broadleaf Forests",
    "Deciduous Broadleaf Forests",
    "Mixed Forests",
    "Closed Shrublands",
    "Open Shrublands",
    "Woody Savannas",
    "Savannas",
    "Grassslands",
    "Permanent Wetlands",
    "Croplands",
    "Urban and Built-up Lands",
    "Cropland/Natural Vegetation Mosaics",
    "Barren",
    "Water Bodies"
  )
)

Please bear in mind that your landcover layer might not have the same classes/levels as mine. Check here the classification values: MCD12_User_Guide_V6 (see page 7 Classification Legends).

Peru: drivers of amphibians richness 🐸

  1. Visualise predictors

Land-cover (MODIS Terra MCD12Q1)

Peru: drivers of amphibians richness 🐸

  1. Visualise predictors

Human footprint index (SEDAC)

Peru: drivers of amphibians richness 🐸

  1. Calculate the correlation between the different predictors (and compare with biodiversity)

Let’s get values per grid-cell.

names(tavg)
 [1] "PER_wc2.1_30s_tavg_1"  "PER_wc2.1_30s_tavg_2"  "PER_wc2.1_30s_tavg_3" 
 [4] "PER_wc2.1_30s_tavg_4"  "PER_wc2.1_30s_tavg_5"  "PER_wc2.1_30s_tavg_6" 
 [7] "PER_wc2.1_30s_tavg_7"  "PER_wc2.1_30s_tavg_8"  "PER_wc2.1_30s_tavg_9" 
[10] "PER_wc2.1_30s_tavg_10" "PER_wc2.1_30s_tavg_11" "PER_wc2.1_30s_tavg_12"

Peru: drivers of amphibians richness 🐸

  1. Calculate the correlation between the different predictors (and compare with biodiversity)

Let’s see all the values together.

predictors <- tibble(
  "tavg" = values(mean(tavg)) %>% as.vector(),
  "prec" = values(mean(prec)) %>% as.vector(),
  "elev" = values(elev) %>% as.vector(),
  "land" = values(land) %>% as.vector(),
  "footp" = values(footp) %>% as.vector()
) %>% 
  mutate(
    across(everything(), ~replace_na(.x, NA))
  )

Peru: drivers of amphibians richness 🐸

  1. Calculate the correlation between the different predictors (and compare with biodiversity)

Let’s see all the values together.

summary(predictors)
      tavg             prec             elev              land       
 Min.   :-0.976   Min.   :  0.00   Min.   : -10.43   Min.   : 1.000  
 1st Qu.:13.301   1st Qu.: 65.48   1st Qu.: 185.44   1st Qu.: 2.000  
 Median :23.433   Median :140.80   Median : 516.08   Median : 2.000  
 Mean   :19.454   Mean   :133.26   Mean   :1491.72   Mean   : 6.152  
 3rd Qu.:25.979   3rd Qu.:197.03   3rd Qu.:2936.59   3rd Qu.:10.000  
 Max.   :27.172   Max.   :490.53   Max.   :5248.80   Max.   :17.000  
 NA's   :16468    NA's   :16468    NA's   :16804     NA's   :16133   
     footp       
 Min.   : 0.000  
 1st Qu.: 1.000  
 Median : 2.000  
 Mean   : 3.316  
 3rd Qu.: 5.000  
 Max.   :43.000  
 NA's   :16133

Peru: drivers of amphibians richness 🐸

  1. Calculate the correlation between the different predictors (and compare with biodiversity)

Let’s see how these values correlate.

GGally::ggpairs(predictors, cex = 0.1)

What can this plot tell you?

EXTRA Can we see any pattern with the hotspots of amphibians species richness?

Peru: drivers of amphibians richness 🐸

Let’s read the data, plot it and see the patterns

amphibians_grid <- readRDS("data/amphibians_grid.rds") %>%
  st_transform(., 4326)

Peru: drivers of amphibians richness 🐸

Let’s read the data, plot it and see the patterns

Peru: drivers of amphibians richness 🐸

Species richness vs mean annual temperature

Peru: drivers of amphibians richness 🐸

Species richness vs precipitation seasonality

Peru: drivers of amphibians richness 🐸

Species richness vs elevation

Peru: drivers of amphibians richness 🐸

Species richness vs land-use

Peru: drivers of amphibians richness 🐸

Species richness vs human footprint

Peru: drivers of amphibians richness 🐸

  1. Get a map of Peru
  2. Get 5 different drivers (raster)
  3. Process all the layers to the same resolution and the extent of Peru
  4. Visualise predictors
  5. Calculate the correlation between the different predictors (and compare with biodiversity)


We are done! No it’s your turn :)

EXCERCISE France: drivers of trees richness

France: drivers of trees richness patterns 🌳

  1. Get a map of France
  2. Get 5 different drivers (raster)
  3. Process all the layers to the same resolution and the extent of France
  4. Visualise predictors
  5. Calculate the correlation between the different predictors (and compare with biodiversity)

France: drivers of trees richness patterns 🌳

This time, you will use the following drivers:

  • Temperature seasonality (WorldClim BIO4)
  • Annual precipitation (WorldClim BIO12)
  • Elevation (SRTM)
  • Land-cover (MODIS Terra MCD12Q1)
  • Human footprint index (SEDAC)

Any doubts?

https://github.com/petrkeil/Spatial-Ecology-and-Macroecology